/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package bd.org.apache.commons.math.distribution;

import bd.org.apache.commons.math.analysis.UnivariateFunction;
import bd.org.apache.commons.math.analysis.solvers.UnivariateSolverUtils;
import bd.org.apache.commons.math.exception.NotStrictlyPositiveException;
import bd.org.apache.commons.math.exception.NumberIsTooLargeException;
import bd.org.apache.commons.math.exception.OutOfRangeException;
import bd.org.apache.commons.math.exception.util.LocalizedFormats;
import bd.org.apache.commons.math.random.RandomDataImpl;
import bd.org.apache.commons.math.util.FastMath;

import java.io.Serializable;

/**
 * Base class for probability distributions on the reals.
 * Default implementations are provided for some of the methods
 * that do not vary from distribution to distribution.
 *
 * @version $Id: AbstractRealDistribution.java 1244107 2012-02-14 16:17:55Z erans $
 * @since 3.0
 */
public abstract class AbstractRealDistribution implements RealDistribution, Serializable {

    /**
     * Default accuracy.
     */
    public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;

    /**
     * Serializable version identifier
     */
    private static final long serialVersionUID = -38038050983108802L;

    /**
     * RandomData instance used to generate samples from the distribution.
     */
    protected final RandomDataImpl randomData = new RandomDataImpl();

    /**
     * Solver absolute accuracy for inverse cumulative computation
     */
    private double solverAbsoluteAccuracy = SOLVER_DEFAULT_ABSOLUTE_ACCURACY;

    /**
     * Default constructor.
     */
    protected AbstractRealDistribution() {}

    /**
     * {@inheritDoc}
     * The default implementation uses the identity
     * <p>{@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}</p>
     */
    @Override
    public double cumulativeProbability(double x0, double x1) throws NumberIsTooLargeException {

        if (x0 > x1) {

            throw new NumberIsTooLargeException(
                LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, x0, x1, true);
        }

        return cumulativeProbability(x1) - cumulativeProbability(x0);
    }


    /**
     * {@inheritDoc}
     * The default implementation returns
     * <ul>
     * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
     * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
     * </ul>
     */
    @Override
    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {

        /*
         * IMPLEMENTATION NOTES
         * --------------------
         * Where applicable, use is made of the one-sided Chebyshev inequality
         * to bracket the root. This inequality states that
         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
         * mu: mean, sig: standard deviation. Equivalently
         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
         * F(mu + k * sig) >= k^2 / (1 + k^2).
         *
         * For k = sqrt(p / (1 - p)), we find
         * F(mu + k * sig) >= p,
         * and (mu + k * sig) is an upper-bound for the root.
         *
         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
         * F(mu - k * sig) <= 1 / (1 + k^2).
         *
         * For k = sqrt((1 - p) / p), we find
         * F(mu - k * sig) <= p,
         * and (mu - k * sig) is a lower-bound for the root.
         *
         * In cases where the Chebyshev inequality does not apply, geometric
         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
         * the root.
         */
        if ((p < 0.0) || (p > 1.0)) {
            throw new OutOfRangeException(p, 0, 1);
        }

        double lowerBound = getSupportLowerBound();

        if (p == 0.0) {
            return lowerBound;
        }

        double upperBound = getSupportUpperBound();

        if (p == 1.0) {
            return upperBound;
        }

        final double  mu  = getNumericalMean();
        final double  sig = FastMath.sqrt(getNumericalVariance());
        final boolean chebyshevApplies;

        chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) || Double.isInfinite(sig)
                             || Double.isNaN(sig));

        if (lowerBound == Double.NEGATIVE_INFINITY) {

            if (chebyshevApplies) {
                lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
            } else {

                lowerBound = -1.0;

                while (cumulativeProbability(lowerBound) >= p) {

                    lowerBound *= 2.0;
                }
            }
        }

        if (upperBound == Double.POSITIVE_INFINITY) {

            if (chebyshevApplies) {
                upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
            } else {

                upperBound = 1.0;

                while (cumulativeProbability(upperBound) < p) {

                    upperBound *= 2.0;
                }
            }
        }

        final UnivariateFunction toSolve = new UnivariateFunction() {

            @Override
            public double value(final double x) {

                return cumulativeProbability(x) - p;
            }
        };

        double x = UnivariateSolverUtils.solve(toSolve, lowerBound, upperBound,
                       getSolverAbsoluteAccuracy());

        if (!isSupportConnected()) {

            /* Test for plateau. */
            final double dx = getSolverAbsoluteAccuracy();

            if (x - dx >= getSupportLowerBound()) {

                double px = cumulativeProbability(x);

                if (cumulativeProbability(x - dx) == px) {

                    upperBound = x;

                    while (upperBound - lowerBound > dx) {

                        final double midPoint = 0.5 * (lowerBound + upperBound);

                        if (cumulativeProbability(midPoint) < px) {
                            lowerBound = midPoint;
                        } else {
                            upperBound = midPoint;
                        }
                    }

                    return upperBound;
                }
            }
        }

        return x;
    }


    /**
     * Returns the solver absolute accuracy for inverse cumulative computation.
     * You can override this method in order to use a Brent solver with an
     * absolute accuracy different from the default.
     *
     * @return the maximum absolute error in inverse cumulative probability estimates
     */
    protected double getSolverAbsoluteAccuracy() {

        return solverAbsoluteAccuracy;
    }


    /**
     * {@inheritDoc}
     */
    @Override
    public void reseedRandomGenerator(long seed) {

        randomData.reSeed(seed);
    }


    /**
     * {@inheritDoc}
     * The default implementation uses the
     * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
     * inversion method.
     * </a>
     */
    @Override
    public double sample() {

        return randomData.nextInversionDeviate(this);
    }


    /**
     * {@inheritDoc}
     * The default implementation generates the sample by calling
     * {@link #sample()} in a loop.
     */
    @Override
    public double[] sample(int sampleSize) {

        if (sampleSize <= 0) {
            throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES, sampleSize);
        }

        double[] out = new double[sampleSize];

        for (int i = 0; i < sampleSize; i++) {

            out[i] = sample();
        }

        return out;
    }
}
